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Abstract 

We calculate the cross section of the exclusive process e + + e~ — > J/ip + g c at the leading order 
approximation within the QCD light-cone sum rules approach. It is found that the form factor 
F\rp{V = J/ip,P = rj c ) depends mainly on the behavior of the twist-2 distribution amplitude of 
the rfc-meson at the scale of this process. Thus in order to obtain a reliable estimation of the cross 
section, it is important to have a realistic distribution amplitude of the rj c meson, and to deal with 
the evolution of the distribution amplitude to the effective energy scale of the process. Our re- 
sults show that one can obtain a compatible prediction with the Belle and BaBar experimental data. 
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I. INTRODUCTION 



The double- charmonium production in e + e~ annihilation at the B factory provides a 
good platform to study both perturbative and non-perturbative effects in quantum chromo- 
dynamics (QCD). 

On one hand, from the experimental point of view, the cross section of the precess e + + 



show that 



point 

Is 



e — > J/ip + rjr was measured at Belle |2] and BaBar [3], and their recent observation 



cx(e + + e~ J/ip + rj c ) x B >2 = 25.6 ± 2.8 ± 3.4 fb (Belle) 



and 



a(e + + e~ ->■ J/ip + rj c ) x B >2 = 17.6 ± 2.8+^ fb (BaBar), 



where B >2 denotes the branching fraction for the final states with more than two charged 
tracks. 

On the other hand, from the theoretical point of view, the process e + + — > J/ip + r] c is 
usually studied within the nonrelativistic QCD (NRQCD) [4]. Under the leading order (LO) 
NRQCD calculation, Refs.{5-7] derived a much smaller cross section a = 3.78 ~ 5.5 fb for 
the first time. Thus there are large discrepancy between theoretical predictions of the LO 
NRQCD calculation and the experimental measurements. In order to solve this problem, 
many attempts were made based on the NRQCD approach. Refs.0, [9] diminished the 
disagreement to a large degree by including the radiative correction. Moreover, as pointed 
out in Refs. 10Ml2l|. the relativistic corrections can further improve the accuracy. By taking 



both the radiative and relativistic corrections into account, Ref.[12( got 17.6lg 7 fb, then the 
authors there optimistically concluded that the disagreement between theoretical predictions 
in NRQCD approach and experiments has been resolved [12]. However, one may doubt the 
validity of the a s - and v 2 - expansion in the process e + + e~ — > J/tp + r] c , if the LO result is 
an order of magnitude smaller than the experimental measurements and the next-to-leading 
order (NLO) corrections/higher t> 2 -expansion terms inversely play a dominate role for the 
double charmonium production. Furthermore, Ref. 9| showed that the scale dependence of 
the cross section can not be improved even with the NLO correction, so it is an important 
matter to determine the typical scale of the process or at least to make a more reliable 
estimation of scale dependence. 
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In contrast, it was argued that the experimental results by Belle and BaBar collaborations 
of the process e + + e~ — > J/ip + r\ c can also be explained by using the perturbative QCD 
(pQCD) with proper models for the charmonium distribution amplitudes (DA) 13) . They 
claimed, "the difficulties in explaining the Belle and BaBar results for cr(e + + e~ — > J/tp + rj c ) 
are not really the difficulties of QCD, but are rather due to a poor approximation of the real 
dynamics of c-quarks by NRQCD". Actually, the exclusive process can be factorized into 
two parts in the pQCD approach: the calculable hard-parton amplitude and the hadronic 
distribution amplitude. If one replaces all DAs in the pQCD formulae by a simple 6 function, 
the calculated cross section will be back to a few fb that is consistent with the NRQCD 
approach. However one always assume that the hadronic distribution amplitude is not non- 
relativistic. Therefore, the key point is that the relativistic DA, instead of the 5 function, 
enhances the cross section of the process e + + e~ — >• J/ijj + r\ c . 

Furthermore, since the QCD light-cone sum rules (LCSR) combines the QCD sum rules 
and the pQCD theory of hard exclusive processes in a suitable way, it can be a good tool 
for calculating the form factors in the large momentum transfer. For example, LCSR is a 
successful method for dealing with the 7*p — > tc transition form factor [14j that is similar 
to the process e + + e~ — > J/ip + r\ c . We shall try to apply the QCD LCSR approach to 
calculate the amplitude of the process e + + e~ — > V + P, where V = J /ip, ... and P = r] c , .... 
The theoretical predictions at the LO approximation for the cross section of the process 
e + + e~ V + P at the large energy scale can be obtained by using the QCD LCSR. Similar 
to the pQCD approach, this method also faces two problems: which charmonium DA model 
should be adopted and how large effects can be determined by the DA evolution. In the 
present paper, we shall discuss the behavior of the different models for the charmonium DA 
and the effects of the renormalized group evolution with the effective scale in the process. 

The remaining parts of the paper are organized as follows. In Sec. II, we present calcula- 
tion technology for e + + e~ — > J/ip + r] c under the QCD LCSR, where the leading-twist DA 
is constructed and its QCD evolution is presented. Numerical results for the cross section 
of the process e + + e~ —> J ' /ip + r] c are presented in Sec. III. The final section is reserved for 
summary and conclusion. 
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II. CALCULATION TECHNOLOGY FOR e+ + e~ -> J/V> + % 



Generally, the cross section for the process a(pi) + 6(^2) — > c(p3) + c?(p 4 ) is given by 
1 f d 3 p 3 G? 3 p 4 



(7 



(27T)V(p 1 +p 2 -p 3 -p4)|Ath 



(1) 



AE 1 E 2 v rd J (2it) 3 2E 3 (2it) 3 2E 4 
where pi = (Ei,pi) (i = 1, • • • ,4) stand for the four-momentum of initial and final particles 
correspondingly, v re i — ||r — |jr|. |.M| 2 is the squared absolute value of the matrix element, 
where the color states and spin projections of the initial and final particles have been summed 
up and those of the initial particles have been averaged. 

For the exclusive double-charmonium production e + +e~ — > J/if)+rj c , its Lorentz-invariant 
matrix element turns out to be 



M=i rf 4 a;<J/^ c |T{Q c j;(x)^(x),e(0) 7 ,e(0)^(0)} 



(2) 



where JZ(x) = C(x)7 M C(x) is the c-quark electromagnetic current. Then, we obtain 



\M\ 



l2 < s - 4m V [lW()] 



VP 



45 



(3) 



where 6 is the scattering angle, Q c = | is the charm quark charge and the form factor Fyp 



is defined as 



<JMP - q) Vc (P)\J c J0)\0> = e mhc e a *q h P c Fv P 



(4) 



with e a being the polarization vector of J/if>- meson and S = —q 2 . Neglecting the small mass 
difference between J /if} and r] c mesons, the cross section becomes 



a 



7ra 2 Q 



2/02 

c I 1 



Am 2 , \ 3/2 

J /i> 1 1 171 12 



6 V " 3 J 

It is shown that the main part is to calculate the form 



VP 



(5) 



methods to calculate it, such as NRQCD [5H12J, PQCD [ljand light -cone perturbative 



QCD approaches [l5|, ll6|]. Here we use the LCSR approach |l7H20j] to calculate Fyp 



actor Fyp- There are many 



A. the form factor F VP within the QCD LCSR 

We adopt the following two-point correlator to calculate the form factor Fyp 
n^(P,g) = i J d 4 xe-^<r ?c (P)|T{^(x)J, c (0)}|0>, 



(6) 



where q is the four-momentum of the virtual photon, P is the four-momentum of f] c meson. 
On one hand, by inserting a complete set of intermediate hadronic states in Eq.(j6]), we 

get 

1 1 f°° ImF 

where the decay constant /j/^ is defined as, <0| J^(0)\ J j^{P — q)> = fj/^mj/^ii, with e M 
being the polarization vector of J/ip- meson. so is the threshold parameter whose value can 
be taken as 3.6 2 GeV 2 < s < 4.2 2 GeV 2 [21]. The second term in Eq.([7j) is the dispersion 
integral that includes the contributions from the excited and continuum states in the region 

S > Sq. 

On the other hand, the correlation function Eq.((6]) can also be calculated by expanding the 
T-product of quark currents near the light cone x 2 = due to sufficiently large momentum 
transfer. For such purpose, we contract the two c-quark fields and write down a free c-quark 
propagator 

CimO) = *(*,0) = i I (8) 
Then up to twist-3 accuracy, Eq.Q can be simplified as 

where m c is the current c-quark mass, f Vc is the decay constant of rj c meson, <f) Vc stands for 
the f] c leading-twist DA that is defined through the matrix element: 



< Vc (P)\C(x)Yl 5 C(0)\V> 



4P T f Vc due mPx (j) ric (u) + higher twist terms. (10) 
Jo 



Next, b y a pplying the quark-hadron duality to Eq.flT]) and by applying the Borel trans- 



formation 



22] 



1 1 -ZlM. 

m - (? - p ) M 
1 1 



M m 2 c -(q-xP) 2 ' xM 2 ' 1 ' 

to Eq.Q and Eq.©, we obtain the sum rule for Fyp 



VP 



2_U /^^J"^^ 1 "^^ 1 ^ 9 ^^] (12) 



m J/-4>fj/^ J A x 



TABLE I: The explicit form of the spin-space wave function x XlXi ( x i k±), where the transverse 
momentum k± = (k x ,k y ) and m* stands for the constituent c-quark mass. 



A1A2 


tt 


tt 


tt 


tt 




kx iky 




m* 


kx~\~iky 


v/2(m^ + fe'|) 


yj2(mf+k'] ) 


V2(m*2+fc2) 





where A 



2i< 



M 1 is the Borel 



{s ~ m\ + Q 2 ) 2 + A(m 2 c + Q 2 )m 2 c - (s - m\ - q 2 
transformation parameter and P 2 = m 2 c , —q 2 = Q 2 = S = 112GeV 2 . It is found that the 
form factor Fyp depends heavily on the DA rf> Vc , especially on its end point behavior due to 



A ~ 0.9. 



B. leading-twist DA of rj c meson 



The key input for the form factor is the gauge-independent and process-independent DA 
which is of non-perturbative nature and can be defined as the integral of the valence 



Fock wave function 



23] 



2^6 



d 2 k ± T , r , 

l07T d 



(13) 



fvc J\k]_\ 2 <^ 

where /io stands for the separation scale between the perturbative and non-perturbative 
regions. As for a scale fi > /Iq, the non-perturbative DA (f) Vc (x,[i) is given by the renormal- 
ization group evolution that can be calculated perturbatively. 

Up to now, it is difficult to give the light-cone wave function (LCWF) from the first 
principles of QCD. So one usually constructs some phenomenological models for the wave 
function such as BC model |r3j, BKL model [2J , BLL model [25j, MS model [2d], BHL 
model 23] and etc. Here, we shall take the BHL model for the rj c wave function, which can 



be written as 



28] 



A1A2 1 

Vc 



x, kj_) = (Pbkl(x, k±)x XlX2 (x, kx) = Ae 



x XlX2 {x, k±), 



(14) 



where m* stands for the constituent c-quark mass, Ai and A2 are helicity states of the 
constitute c and c quarks, x XlX2 { x i k±) stands for the spin-space wave function coming from 



the Wigner-Melosh rotation 29|. x AlA2 (x, kj_) can be found in Refs. 30M32[| . whose explicit 
form is shown in TABlU 
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FIG. 1: The comparison of the i] c DAs of BHL model with those of BC [13j, BKL (240, BLL [25] 
and BHL models at the initial scale fj,Q. 

The parameters A and b 2 can be determined by two constraints. One constraint is from 
the wave function normalization 

^/W ^£n iA2 (^) = l- (15) 
U Jo J\kl\'<fi 16tt 3 ^ 

It can be found that only the usual helicity components Ai + A2 = makes contribution to 
the wave function normalization. More explicitly, from TAB]U we have x Al+A2=0 (x, k±) = 
, ™ a :. Another constraint is from the probability of finding the leading Fock state |cc > 
in the r\ c Fock state expansion, i.e. 

1 r d 2 k 

dx I y^|y?BHL(>, k ± )\ 2 = P Vc , (16) 







withP„ ~0.8 

'/c 



30|. One can assume [Xq = m* to be the initial scale for the non-perturbative 



distribution amplitude of the rj c - meson. Inputting the constituent quark mass m* = 1.8 GeV 



[131 ]. the decay constant f Vc = 0.335 GeV [33J and the initial scale fi = m* = 1.8 GeV, we get 
the corresponding parameters for A = 285.64291 GeV _1 ,6 2 = 0.19057 GeV -2 . We compare 



the r? c -DA of our BHL model at the scale //„ with those of BC Q, BKL |24j, BLL [25] 
models in Fig.(pQ). 

However, the scale /1 of the form factor (|T2|) that is at the B-factory is very different from 
the initial scale /j,q of rj c DA. When DA runs to a higher energy scale fj,, other than fi , with 
a proper QCD evolution, the behavior of DA shall be changed to a certain degree, especially 
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in its upper end-point regions A < x < 1 that determines the form factor Fyp as shown by 
Eq. ffT2l) . Therefore, it is quite important to do the DA evolution from the initial scale /io 
to a typical energy scale /i of the process so as to derive a more reliable cross section for 
e + + e~ — > J /if) + rj c process. Thus, the next section is devoted to deal with the evolution of 
the i] c DA. 



C. The evolution of the y c DA with the scale fi 



We describe the DA evolution according to Ref . 23] . In the light-cone gauge, the DA (f) nc 
is related to the hadronic wave function ^/ Vc , which is the Fourier transform of the positive- 
energy projection of the usual Bethe-Salpeter wave function evaluated at relative "light-cone 
time", i.e. 

,2 \ -7f//9 



1 ^ 

ln To 
A 2 



d 2 k ± 

16vr 3 



(17) 



where the factor in front of the integral comes from the scale dependence due to vertex and 
self-energy insertions. An evolution equation is obtained by differentiating both sides of 
Eq. ljTTI) with respect to jj 2 . To order 0(a s ), we obtain an "evolution equation" [23] 



XlX 2 /J, = = C F — - — 



[dy]V(x i ,y l )(p Vc (yi,fi) - xix 2 4> Vc (xi, lA 



18) 



where 



xiy 2 0(yi - xi) 5 hlfl2 + 



V{x hyi ) = 2C F 

[dy] = dy 1 dy 2 5(l -yi- y 2 ) 

4>Voi X ii A*) = XlX2<t>ri c {Xi,ti), 



A 



(yi - xi) 



+ (1 O 2) 



(19) 



Cp = 4/3, 8 hl h 2 = 1 when the c and c helicities are opposite, and A0 J?c (?/ i , /i) = 4> Vc (yi, fi>) 
(firjcixi, /i). The running coupling constant a s (/x 2 ) at the LO is given by a s (/i 2 ] 



&o in 

with &o — 25/3. One explicit solution of Eq. (|T8|) can be written in the following Gegenbauer 
expansion 

,2 \ ~ln 



n=0 



, c (xi,/i) = x x x 2 ^a n ( ln^ ) C^ 2 {xi-x 2 ) 



(20) 
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FIG. 2: rj c -T)A derived with strict evolution (118p at some typical energy scales, where the left is for 



BHL model 



3 



and the right is for BLL model 25j. The solid lines, dashed lines and the dotted 



lines represent DA at fj, = /xo, M = 3.46 GeV and [i = 5.00 GeV respectively. 

where the Gegenbauer polynomials Cn 2 are eigenfunctions of V(xi, yi) and the corresponding 
eigenvalues are the "non-singlet" anomalous dimensions 

*-T( 1 + 4 S£-5H^)* a (21) 

The coefficients a n which are non-perturbative can be determined from the initial condition 
4> Vc (xi, fi ) by using the orthogonality relations for the Gegenbauer polynomials Cn 2 . 

Usually, one truncates the Gegenbauer expansion (120 j) with the first 3 or 4 terms (n = 
0,2,4,6 in our case) to obtain the behavior of DA at the higher energy scales. In this 
paper we solve the evolution equation Eq.f lTB"]) strictly to get the DA's behavior at the large 
scale since DA's behavior is very important for calculating the form factor of the process 
e + + e~ — > J/ip + r] c . The Eq. ffTB"]) and Eq. fl20|) are equivalent to each other if the Gegenbauer 
expansion converges quickly. The evolution of DA with the strict evolution fll8p are shown 
in Fig. ([5]), where the solid lines represent the DAs at the initial energy scale /io = 1.80 GeV, 
the dashed lines and dotted lines represent the DAs at the energy scale = 3.46 GeV and 



/i = 5 GeV that are taken by Ref.[13| and Ref.|25[ respectively. It is shown that when 
the energy scale becomes larger, the DA becomes lower in the middle while becomes higher 
near the end point, till at last when the energy scale tends to infinity, the DA tends to a 
asymptotic form 4> as (x) = 6x(l — x). 
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FIG. 3: The dependance of cross section on the threshold parameter so within the LCSR approach. 
The left and the right correspond to the scale /i = 3.46 GeV and fi = 5.00 GeV with the BHL 
model of rj c DA. 

III. NUMERICAL RESULTS AND DISCUSSION 

To calculate the form factor and the cross section of e + + e~ — > J/ip + r) c , we take 



fj/^p = 0.416 GeV, f Vc = 0.335 GeV, m J/f = 3.096916 GeV and m Vc = 2.9798 GeV [34j,[35|. 
And to compare with the results in literature we take the c-quark current mass to be 
m c = 1.2 GeV. The Borel parameter M 2 ranges from 8 to 15 GeV 2 , when at this range, 
both the form factor and the cross section are stable. As for the effective scale /i of the 
process e + + e~ — > J/ip + 7] c , Ref.[3] suggested jj, ~ Vk 2 w 3.46 GeV from the mean value of 
< Z 1 ^ 0.80 or from the coupling constant < a s (k 2 ) >^ 0.263. Another usually adopted 
scale is » y/S/2 ~ 5 GeV [25J. Here, we will take fi = 3.46 GeV and /i = 5.00 GeV to do 
our discussion. 

As for the threshold parameter s , Ref.[2l] took 3.6 2 GeV 2 < s < 4.2 2 GeV 2 [21] 



with the central value so = 3.8 2 GeV as their case. Similarly, we also take so within the 
same region while with a little different central value. To see the dependence of the cross 
section on the threshold parameter so, we give the cross section corresponding to so = 
3.7 2 , 3.9 2 , 4.1 2 GeV 2 with the BHL model at the scale n = 3.46 GeV and /i = 5.00 GeV in 
Fig-©. Since the cross section with the threshold parameter sq = 3.9 2 GeV 2 is more stable 
than that with sq = 3.7 2 GeV 2 , 4.1 2 GeV 2 , we take sq = 3.9 2 GeV 2 as our central value of 
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FIG. 4: The cross section of e + + e~ — > J/ip + r\ c at different effective scale within the LCSR 
approach. The solid lines, dashed lines and dotted lines correspond to the rj c meson distribution 
amplitude at the scale u = uq, a = 3.46 GeV and u = 5.00 GeV, where the left is for BHL model 

n 



[271 ] and the right is for BLL model 



25]. 



TABLE II: Comparison of the cross section of e + + e 
within the QCD light-cone sum rules approach. 



J/tp + rj c from different r\ c DA models 






BHL 


BHL 


BLL 


BLL 


H (GeV) 


3.46 


5.00 


3.46 


5.00 


a(fb) 


13.08± 0.32 


25.96±0.55 


10.34±0.17 


21.16±0.34 



threshold parameter. 

It is found that the cross section depends on the DA in the region x > A through the 
form factor formula (jl~2l) . We show the cross section corresponding to three typical scales 
(j, = ^ = 1.80 GeV, fj, = 3.46 GeV and [i = 5.00 GeV in Fig.®. When the effective 
energy scale increases, the corresponding cross section becomes bigger and is compatible 
with the BaBar and Belle's measurements. Thus, by setting the effective energy scale and 
dealing with the DA evolution properly, the LCSR can provide a possible explanation for the 
double-charmonium production process e + + e~ — > J/ip + r] c at the B factory. To explicate 
the cross sections of different models numerically, we further show these cross sections in 
Tab. [Til where the error is caused by the variation of M 2 . One may observe that the cross 
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section by BLL model is smaller than that by BHL model since the DA of the BLL model 
is narrower than that of the BHL model as shown in Fig. ([2]). 

The above calculation is done at the LO approximation within the QCD LCSR and the 
error is only caused by the variation of the Borel parameter M 2 . Of course, one should 
include the higher order contributions, such as the NLO corrections to the light-cone sum 
rules, the higher-twist DAs, the higher Fock states and etc. Therefore we can not estimate 
all of the uncertainties of the calculated cross section before doing a further study. 

IV. SUMMARY 

The exclusive charmonium J/ip + i] c production in e + e~ collision is a very interesting 
problem. Since the discrepancy between theoretical prediction of the LO NRQCD and 
experimental data given by Belle and BaBar at the B factory has posed a significant chal- 
lenge for several years, many theoretical attempts have been made to solve this challenging 
problem. It is worthwhile to study this process by taking various applicable approaches to 
understand the charmonium production dynamics. In this paper we study this process by 
using the QCD LCSR approach. 

Our results based on the LCSR approach shows that the cross section of the process 
e + + e~ — > J/ip + r] c substantially depends on the behavior of the rj c DA at the energy scale 
/i. Noticing that the energy scale /i at the B factory is greater than the initial scale /io of the 
r] c DA, the renormalization group evolution of the DA has to be taken into account. The 
perturbative radiative correction leads to a big change of the DA especially to the tail of 
the DA at the large scale \i. 

At the present, one has poor knowledge of the DA and tries to build various models that 
have quite different behavior, especially at the end-point region. We stress that the evolution 
of the rj c DA can give more reasonable prediction to the process e + + e~ — > J/ip + rj c within 
the LCSR approach. Similar to other approaches, in order to calculate the cross section of 
the double- charmonium production one needs to have more knowledge of the charmonium 
DA. 

The numerical results show that the cross section of the process e + + e~ — > J/ip + r] c 
is predicted in the range 13 ~ 26/6. The calculated values for the different models can be 
compatible with the Belle and BaBar measurements by properly choosing effective energy 
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scale for this process and dealing with the DA evolution effect. 
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